import java.util.*;
import java.io.*;

public class KernelDensityEstimator {

    public KernelDensityEstimator(final double[] x, final double h, final int n, final double[] lims, boolean bandwdithLimited) {
	this.bivariate = false;
        this.x = x;
	this.y = x;
        this.nx = x.length;
        if (n <= 0) throw new RuntimeException("must have a positive number of grid points");
        this.n = n;
        if (lims != null) this.lims = lims; else setupLims();
        this.limitBandwidth = bandwdithLimited;
        if (h >= 0) this.h = new double[]{h}; else setupH();
        doKDE1D();
    }

    public KernelDensityEstimator(final double[] x, final double h, final int n, final double[] lims) {
        this(x, h, n, lims, false);
    }

    public KernelDensityEstimator(final double[] x) {
        this(x,-1.0,50,null);
    }

    public KernelDensityEstimator(final double[] x, final int n) {
        this(x,-1.0,n,null);
    }

    public KernelDensityEstimator(final double[] x, final double[] y, final double[] h, final int n, final double[] lims) {
        this(x, y, h, n, lims, false);
    }

    public KernelDensityEstimator(final double[] x, final double[] y, final double[] h, final int n, final double[] lims, boolean bandwdithLimited) {
	this.bivariate = true;
        this.x = x;
        this.y = y;
        if (x.length != y.length) throw new RuntimeException("data vectors must be the same length");
        this.nx = x.length;
        if (n <= 0) throw new RuntimeException("must have a positive number of grid points");
        this.n = n;
        if (lims != null) this.lims = lims; else setupLims();
        this.limitBandwidth = bandwdithLimited;
        if (h != null) this.h = h; else setupH();
        doKDE2D();
    }
    
    public KernelDensityEstimator(final double[] x, final double[] y, boolean limitBandwidth) {
        this(x,y,null,50,null,limitBandwidth);
    }

    public KernelDensityEstimator(final double[] x, final double[] y) {
        this(x,y,null,50,null);
    }

    public KernelDensityEstimator(final double[] x, final double[] y, final int n) {
        this(x,y,null,n,null);
    }

    public void doKDE2D() {
	if ( bivariate == false ) throw new RuntimeException();
        gx = makeSequence(lims[0], lims[1], n);
        gy = makeSequence(lims[2], lims[3], n);
        double[][] ax = outerMinusScaled(gx, x, h[0]);
        double[][] ay = outerMinusScaled(gy, y, h[1]);
        normalize(ax);
        normalize(ay);
        z2D = new double[n][n];
        double scale = nx * h[0] * h[1];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                double value = 0;
                for (int k = 0; k < nx; k++) value += ax[i][k] * ay[j][k];
                z2D[i][j] = value / scale;
            }
        }
    }

    public void doKDE1D() {
	if ( bivariate == true ) throw new RuntimeException();
        gx = makeSequence(lims[0], lims[1], n);
        double[][] ax = outerMinusScaled(gx, x, h[0]);
        normalize(ax);
        z1D = new double[n];
        double scale = nx * h[0];
        for (int i = 0; i < n; i++) {
                double value = 0;
                for (int k = 0; k < nx; k++) value += ax[i][k];
                z1D[i] = value / scale;
        }
    }

    public double findLevelCorrespondingToMass2D(double probabilityMass) {
	if ( bivariate == false ) throw new RuntimeException();
        double level = 0;
        double[] sz = new double[n*n];
        double[] c1 = new double[n*n];
        for(int i=0; i<n; i++) System.arraycopy(z2D[i],0,sz,i*n,n);
        Arrays.sort(sz);
        final double dx = gx[1] - gx[0];
        final double dy = gy[1] - gy[0];
        final double dxdy = dx * dy;
        c1[0] = sz[0] * dxdy;
        final double criticalValue = 1.0 - probabilityMass;
        if (criticalValue < c1[0] || criticalValue >= 1.0) throw new RuntimeException();
        for(int i=1; i<n*n; i++) {
            c1[i] = sz[i] * dxdy + c1[i-1];
            if (c1[i] > criticalValue) {
                final double diffC1 = c1[i] - c1[i-1];
                final double diffSz = sz[i] - sz[i-1];
                level = sz[i] - (c1[i]-criticalValue) / diffC1 * diffSz;
                break;
            }
        }
        return level;
    }

    public double[][] getKDE2D() { 
	if ( bivariate == false ) throw new RuntimeException();
	return z2D; 
    }

    public double[] getKDE1D() { 
	if ( bivariate == true ) throw new RuntimeException();
	return z1D;
    }

    public double[] getGrid() { 
	if ( bivariate == true ) throw new RuntimeException();
	return gx; 
    }

    public double[] getXGrid() { return gx; }

    public double[] getYGrid() { 
	if ( bivariate == false ) throw new RuntimeException();
	return gy; 
    }

    public void normalize(double[][] X) {
        for (int i = 0; i < X.length; i++) {
            for (int j = 0; j < X[0].length; j++)
                X[i][j] = normalDistributionPDF( X[i][j], 0, 1);
        }
    }

    public void normalize(double[] X) {
        for (int i = 0; i < X.length; i++) X[i] = normalDistributionPDF( X[i], 0, 1);
    }

    public static double normalDistributionPDF(double x) {
        return Math.exp(-x*x / 2) / Math.sqrt(2 * Math.PI);
    }

    public static double normalDistributionPDF(double x, double mu, double sigma) {
        return normalDistributionPDF((x - mu) / sigma) / sigma;
    }

    public static double normalDistributionCDF(double z) {
        if (z < -8.0) return 0.0;
        if (z >  8.0) return 1.0;
        double sum = 0.0, term = z;
        for (int i = 3; sum + term != sum; i += 2) {
            sum  = sum + term;
            term = term * z * z / i;
        }
        return 0.5 + sum * normalDistributionPDF(z);
    }

    public static double normalDistributionCDF(double z, double mu, double sigma) {
        return normalDistributionCDF((z - mu) / sigma);
    } 

    public double[][] outerMinusScaled(double[] X, double[] Y, double scale) {
        double[][] A = new double[X.length][Y.length];
        for (int indexX = 0; indexX < X.length; indexX++) {
            for (int indexY = 0; indexY < Y.length; indexY++)
                A[indexX][indexY] = (X[indexX] - Y[indexY]) / scale;
        }
        return A;
    }

    public double[] makeSequence(double start, double end, int length) {
        double[] seq = new double[length];
        double by = (end - start) / (length - 1);
        double value = start;
        for (int i = 0; i < length; i++, value += by) seq[i] = value;
        return seq;
    }

    private double margin = 0.1;

    private static double min ( double[] vals ) {
	double min = Double.MAX_VALUE;
	for ( double v : vals ) if ( v < min ) min = v;
	return min;
    }

    private static double max ( double[] vals ) {
	double max = Double.MIN_VALUE;
	for ( double v : vals ) if ( v > max ) max = v;
	return max;
    }

    public static double mean ( double[] data  ) {
        double sum = 0.0;
        for(double a : data) sum += a;
        return sum / ((double)data.length);
    }

    public static double variance ( double[] data ) {
            double mean = mean(data);
            double temp = 0;
            for(double a :data) temp += (mean-a)*(mean-a);
            return temp / ((double)data.length);
    }

    public static double stddev ( double[] data ) {
            return Math.sqrt(variance(data));
    }

    public static double median ( double[] data ) {
               double[] b = new double[data.length];
               System.arraycopy(data, 0, b, 0, b.length);
               Arrays.sort(b);
               if (data.length % 2 == 0) return (b[(b.length / 2) - 1] + b[b.length / 2]) / 2.0;
               else return b[b.length / 2];
    }

    private void setupLims() {
	if ( bivariate ) { 
	        lims = new double[4];
	        lims[0] = min(x);
	        lims[1] = max(x);
	        lims[2] = min(y);
	        lims[3] = max(y);
	        double xDelta = (lims[1] - lims[0]) * margin;
	        double yDelta = (lims[3] - lims[2]) * margin;
	        lims[0] -= xDelta;
	        lims[1] += xDelta;
	        lims[2] -= yDelta;
	        lims[3] += yDelta;
        } else {
	        lims = new double[2];
	        lims[0] = min(x);
	        lims[1] = max(x);
	        double xDelta = (lims[1] - lims[0]) * margin;
	        lims[0] -= xDelta;
	        lims[1] += xDelta;
	}
    }

    private void setupH() {
	if ( bivariate ) { 
	        h = new double[2];
	        h[0] = bandwidthNRD(x) / 4;
	        h[1] = bandwidthNRD(y) / 4;
	        if (limitBandwidth) {
	            if (h[0] >  0.5) h[0] = 0.5;
	            if (h[1] > 0.5)  h[1] = 0.5;
	        }
	} else {
	        h = new double[1];
	        h[0] = 1.06 * stddev(x) * Math.pow(x.length , -0.2 );
	        if (limitBandwidth) {
	            if (h[0] >  0.5) h[0] = 0.5;
	        }
	}
    }

    public static double quantile(final double[] values, final double p) { 
	return quantile(values, 0, values.length, p);
    }

    public static double quantile(final double[] values, final int begin, final int length, final double p) {
        if ((p > 100) || (p <= 0)) throw new IllegalArgumentException("invalid quantile value: " + p);
        double n = (double) length;
        if (n == 0) return Double.NaN;
        if (n == 1) return values[begin];
        double pos = p * (n + 1) / 100;
        double fpos = Math.floor(pos);
        int intPos = (int) fpos;
        double dif = pos - fpos;
        double[] sorted = new double[length];
        System.arraycopy(values, begin, sorted, 0, length);
        Arrays.sort(sorted);
        if (pos < 1) return sorted[0];
        if (pos >= n) return sorted[length - 1];
        double lower = sorted[intPos - 1];
        double upper = sorted[intPos];
        return lower + dif * (upper - lower);
    }

    public double bandwidthNRD(double[] in) {
        final double h = (quantile(in, 75) - quantile(in, 25)) / 1.34;
        return 4 * 1.06 * Math.min(Math.sqrt(variance(in)), h) * Math.pow(in.length, -0.2);
    }

    public double[] getLims() { return lims; }

    public static void main(String[] arg) {
        double[] x = {3.4, 1.2, 5.6, 2.2, 3.1};
        double[] y = {1.0, 2.0, 1.0, 2.0, 1.0};
        KernelDensityEstimator kde = new KernelDensityEstimator(x, y, 4);
        double gridX[]   = kde.getXGrid();
        double gridY[]   = kde.getYGrid();
        double kde2d[][] = kde.getKDE2D();
        KernelDensityEstimator kde2 = new KernelDensityEstimator(x);
        double grid[]  = kde2.getGrid();
        double kde1d[] = kde2.getKDE1D();
    }

    private final double[] x; // x coordinates
    private final double[] y; // y coordinates
    private double[] h;       // h[0] x-bandwidth, h[1] y-bandwidth
    private final int n;      // grid size
    private double[] lims;    // x,y limits
    private int nx;           // length of vectors
    private double[] gx;      // x-grid points
    private double[] gy;      // y-grid points
    private double[][] z2D;   // 2D KDE estimate;
    private double[] z1D;     // 1D KDE estimate;

    private final boolean limitBandwidth;

    private final boolean bivariate; 

}
